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Abstract 

We give an overview of track fitting using the Kalman filter method in the NOMAD 
detector at CERN, and emphasize how the wealth of by-product information can 
be used to analyze track breakpoints (discontinuities in track parameters caused by 
scattering, decay, etc.). After reviewing how this information has been previously 
exploited by others, we describe extensions which add power to breakpoint detection 
and characterization. We show how complete fits to the entire track, with breakpoint 
parameters added, can be easily obtained from the information from unbroken fits. 
Tests inspired by the Fisher F-test can then be used to judge breakpoints. Signed 
quantities (such as change in momentum at the breakpoint) can supplement un- 
signed quantities such as the various chisquares. We illustrate the method with 
electrons from real data, and with Monte Carlo simulations of pion decays. 
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1 Introduction 



The Kalman filter is an efficient algorithm for fitting tracks in particle spec- 
trometers with many position-sensing detectors [1-6]. It cures many of the 
problems of traditional \ 2 track fitting using Newton steps, which becomes 

1 Email address: cousins@physics.ucla.edu 

2 on leave from the Laboratory of Nuclear Problems, JINR, 141980 Dubna, Russia. 



Accepted for publication in Nuclear Instruments and Methods 14 December 1999 



more and more unwieldy as the number of position measurements increases. In 
such situations where Kalman filtering is naturally applied, it can be possible 
to detect and characterize track breakpoints, defined as locations where one or 
more of the track parameters is discontinuous. Obvious breakpoints, such as 
a large kink due to a particle decay, are often found before a full track fit is 
performed. More subtle breakpoints may only manifest themselves when the 
track is fit to obtain the track parameters. 

Fruhwirth[3] has investigated the detection of breakpoints using information 
which is a natural by-product of a Kalman filter track fit. When fitting a drift- 
chamber track with N position measurements ("hits"), the idea has been the 
following. At the location of hit k, one has the best-fit track parameters (a) 
using hits 1 through k and (b) using hits k + 1 through N. One constructs a x 2 
for the consistency of these two sets of track parameters. This x 2 5 along with 
other x 2 's at hand from the two fits, can be combined to form test statistics 
for breakpoints. 

Such methods suffer a loss of power because of two defects. First, x 2 's by con- 
struction throw away information about the arithmetic signs of differences; 
such information is relevant since a track's momentum should normally de- 
crease when the particle decays. 

Second, appropriate constraints of a breakpoint hypothesis are not incorpo- 
rated. For example, if a particle decays at hit k, a desired quantity is the 
mismatch in track parameters describing the momentum vector, under the 
constraint that the track position vector from fits (a) and (b) is identical. 

In this paper, we describe a procedure which uses information from the Kalman 
filter fit to construct the result of a full track fit which has additional track 
parameters to account for the discontinuities at a breakpoint. It is natural 
to allow for one, two, or three discontinuous parameters in order to describe 
different physical processes: 

Type I: An electron emitting a bremsstrahlung photon generally changes 
only its momentum magnitude, since the photon is essentially collinear with 
the electron direction. 

Type II: A particle with a hard elastic scatter may have momentum mag- 
nitude essentially unchanged, while changing the two angles specifying di- 
rection. 

Type III: A charged pion or kaon decaying to /iu in general changes mo- 
mentum magnitude as well as the two angles. 

With the method we describe, trial breakpoint fits at every hit (away from 
the ends) of the track can be quickly obtained, and used to search for and 
characterize breakpoints. 
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Billoir [1] has investigated Type II breakpoints, performing fits which do not 
assume an existing breakpointless fit. We show in this paper how the by- 
products of a Kalman filter fit allow one to avoid refitting the hits while 
incorporating the constraints. 

We discuss these tools in the context of the NOMAD [7] neutrino experiment 
at CERN, within which this development took place. However, the results 
are generally applicable to any experiment in which the number of position 
measurements is large enough to fit the track on both sides of the potential 
breakpoint. 

In Sec. 2 we introduce some notation and describe our track parameters and 
track models, including energy loss. In Sec. 3 a review of the traditional (non- 
Kalman) track fit is given. Sec. 4 describes its replacement by the Kalman 
filter. In Sec. 5, we briefly review previous work on breakpoint variables. In 
Sec. 6, we introduce the new breakpoint variables, and in Sec. 7, we present 
some indicative results of their use. We conclude in Sec. 8. 



2 Track parameters and track models 

2. 1 Parameters 

When fitting a track, one typically describes its location in 6D phase space 
by choosing a fixed reference surface (a vacuum window, chamber plane, etc.) 
and then fitting for the 5 independent parameters of the track at the position 
where the track intersects this surface. We let the vector x contain these track 
parameters. In a fixed-target experiment with beam direction along the z axis, 
the parameterization of x is often taken to be 

x = (x,y, dx/dz, dy/dz } q/p) (1) 

at a reference plane at fixed z, where q/p is the charge/momentum. 

NOMAD is a fixed-target experiment with drift chamber planes perpendicular 
to the z direction (nearly aligned with the neutrino beam). However, the 
chambers are immersed in a uniform magnetic fielaQ so that soft tracks often 
loop back, and a helical parameterization similar to collider experiments is 



3 Sense wires of one chamber make angles of +5, and -5 degrees with respect 
to the magnetic field direction providing a space measurement along coordinates 
designated u, y, and v. 
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Fig. 1. Definition of the helix parameters used to describe charged particle trajectory 
in the NOMAD setup. 

more appropriate. We maintain the reference surface as a plane with fixed z 
(the z of the first hit of the track), and specify a track there by 

x= (x,y, l/ff,tanA,0,t), (2) 



where we have introduced the three parameters of a helical curve in the uni- 
form magnetic field: the signed^] inverse radius of curvature 1/R, the dip angle 
tan A, and the angle of rotation </> (see Fig. 1). In addition to these 5 traditional 
parameters, NOMAD has a sixth parameter, the zero-time-offset for the drift 
chamber measurements, called t. It has been introduced because of trigger 
time jitters. In the following, unless specified otherwise, track parameters will 
refer to this second parameterization. 

The parameter 1/R is related to the momentum p by 

i cc -yl + ttanA) 2 , (3) 



while the sign of {1/R) only reflects the particle charge if the time runs the 
right way along the track. In NOMAD the sign convention we used implied 
that the product R ■ increases with time along the assumed time direction 
given by the ordering of the measurements along the track. 

4 In NOMAD, 1/R has a sign opposite to the particle charge. 
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Throughout our work, we use a change in 1/R as an indicator of a change 
in p; this is strictly true only when the change in tan A is negligible, but is 
an adequate approximation. (Inhomogeneity in the magnetic field B does not 
matter, since we compare 1/R estimates at the same point.) 

With this parameterization, the three physical processes enumerated in the 
Introduction have the following breakpoint signatures: 

Type I: An electron emitting a bremsstrahlung photon has a discontinuity 
in 1 /R. 

Type II: A pion with a hard elastic scatter has discontinuities in tan A and 
0. 

Type III: A charged pion or kaon decaying to jj,u has discontinuities in 1/R, 
tan A, and 0. 

From the estimation of the track parameters at the reference plane one needs 
a transformation, called the track model, with which one computes the ex- 
pected measurements at any position in the detector. This model describes 
the dependence of the measurements on the initial values in the ideal case of 
no measurement errors and of deterministic interactions of the particle with 
matter. Use of a correct description is of the utmost importance for the per- 
formance of the fitting procedure, be it traditional or not. 



2.2 Equation of motion in a magnetic field 



The trajectory of a charged particle in a (static) magnetic field is determined 
by the following equation of motion: 

d 2 r/ d 2 s = (kq/p) ■ ( dr/ ds) x B(r(s)), (4) 



where r is the position vector, s is the path length, A; is a constant of propor- 
tionality, q is the (signed) charge of the particle, p is the absolute value of its 
momentum, and B(r) is the static magnetic field. 

With our parameterization it proved convenient to use as the running param- 
eter rather than z since particles may loop back in the detector and cross the 
same measurement plane several times. The equation of motion can be readily 
integrated along a trajectory step (from position to position 1) where R and 
tan A are assumed constant. 



xi = x + R ■ tan A • (0i — O ) (5) 
V\ =Vo + Ro- (cos 0i - COS 0o ) (6) 
zi = z + R - (sin 0i - sin O ) (7) 
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Ri = R (8) 

tan Ai = tan A (9) 

h =t + R -((f) 1 - ( j) )/((3cosXo). (10) 

In the last equation, (3 is the particle velocity and R ■ A(f> has the right sign 
following our convention. 

In NOMAD, detector measurement planes are located at fixed z and Eqn. 7 
can be solved, 

sin0! = sin0 o + Oi - Zq) • (V^o), (H) 



to obtain <pi at the desired z. In practice, among all the possible solutions, 
our track model returns the one corresponding to the next crossing in the 
requested time direction. 

The magnetic field strength varies by a few percent in the tracking volume. 
This was accommodated by ignoring the minor components of the field and 
updating 1/R at every tracking step so that the product R-B remains constant, 
up to energy losses which are now discussed. 



2.3 Energy losses 



The ionization losses are accounted for by updating 1/R at every tracking 
stepQ: 

A(1/R) = iQM^A, = (12) 
v ' ' d£ dx 0.35/? dx y K ' 



where A0 = <fi en d — <j) s tart, and where dE/ dx is given by the Bethe-Bloch 
equation, evaluated (by default) with the pion mass, and a local matter density 
extracted from the detector model used in the GEANT [8] simulation of the 
experiment. In the central tracker part of the detector, the matter density is 
about 0.1 g/cm 3 : the ionization loss model does not need to include detailed 
relativistic corrections. 

In NOMAD, bremsstrahlung losses should be accounted for (on average) in the 
track model for electron tracks, because the central tracker amounts to about 
1 radiation length (Xq). In the electron (and positron) track model, one adds 

5 The absolute value derives from our sign conventions, so that a track gains energy 
if tracked backwards in time. 
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to the ionization losses (evaluated with the electron mass) bremsstrahlung 
losses: 



A(l/R) 



A<f) 



(13) 



(3 2 X cos A ' 



where we can readily approximate (3 — 1. As expressed, these losses include 
the whole radiated photon spectrum, although beyond a certain threshold, the 
radiated photons can be detected in the downstream electromagnetic calorime- 
ter, or sometimes as a conversion pair in the tracker. But no convincing way 
was found to define a threshold that separates continuous small losses from 
accidental big ones. 



3 Traditional track fits 

Given the track model and an estimate (xo) of the parameters at the refer- 
ence position, one can compute at each measurement location in the detector 
(labeled k = 1, . . . N) the theoretical "ideal measurements" which would be 
made in the absence of fluctuations due to the two categories of "noise" : pro- 
cess noise (multiple scattering, bremsstrahlung, etc.) and measurement noise 
(detector resolution). This computation is represented by the system equation 
in the absence of noise, 



where f is a deterministic function (the track model) giving the values of the 
parameters at each location k. 

In general the set of parameters x^ is not measured directly by the apparatus; 
only a function of it, /ifc(x^), is observed. (In our case is a drift chamber 
position measurement.) Let 



be the measurement equation where represents the measurement errors. 
By convention, (m^) = fofc(xfc) and (£j) = 0. In practice, one has to check 
that the track model (involved in the calculation of X&) and the measurement 
function hk fulfill this convention. 



Xfc = ffc(xo), 



(14) 




(15) 
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In a traditional (non-Kalman) track fit, one calculates the NxN covariance 
matrix of the measurements^: 

v^f = (K - ^(x 4 ))K- - ^( Xi ))>, (16) 

where the angle-brackets represent an average over an ensemble of tracks with 
the same track parameters x. According to the Gauss-Markov theorem, the 
minimization of 

X t 2 rk (x) = (m - h(x)) T [V^]- 1 (m - h(x)) (17) 

yields parameter estimates with minimum variance among all linear unbiased 
estimates. 

Note that when the parameter evolution is not affected by any stochastic noise 
(as assumed by Eq. 14), = (&* Ej) is diagonal (at least block diagonal, if 
the apparatus provides multidimensional correlated measurements), and may 
depend weakly on x via the detector response (for example, the spatial reso- 
lution of drift chambers depends on the track angle w.r.t the anode plane). 

The system equation becomes non deterministic when the track experiences 
stochastic processes such as multiple scattering, bremsstrahlung or ionization 
losses. The system equation (Eq. 14) becomes : 

x fc = ffe(xo) + w fc (18) 

where is a random vector representing the fluctuation of the parameters 
along the path from location to location k. This process noise translates, via 
the track model, into off-diagonal elements in y( m ) (the noise from location 
to k depends on the noise from location to k — 1 and k — 1 to k). 

More importantly, now depends on the reference location z at which the 
parameters will be estimated. Of course the average value of the perturbing 
effects have to be included in but the fluctuations (e.g. track scattering or 
energy loss straggling when relevant) have to be described by the p.d.f of w^. 
In fact only the covariance matrix of the is needed in practice. (See Sec. 4.) 

Given particular data m, the track fit consists in finding the value of x which 
minimizes Eq. 17. Minimization is typically an iterative process with some 
convergence criteria to decide when to stop iterating. The matrix can 
be calculated by Monte Carlo techniques or sometimes analytically. If one is 

6 The superscript in parentheses in ~l/( m ) is to make clear which vector it corre- 
sponds to; this will be the convention in this paper. 
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fortunate, the calculation depends only weakly on x, so it can be done once 
per fit using the initial guess for x, and not changed at each iteration. Note 
that the NxN matrix must be inverted, where N is the number of 

measurement positions. 

Let 

Xtrk (written without explicit argument x) be the minimum value after 
convergence of the minimization procedure and x the value of the estimated 
parameters (the value of x at the minimum). The covariance matrix of the 
estimates is then approximated by the inverse of the curvature matrix at the 
minimum : 



V® = 



1 d\j k 

2 dxjdxi 



(19) 



A traditional fit gives the track parameters x only at the fixed reference zq, say 
at the beginning of the track. In order to find the best-fit track parameters at 
the end of the track, one has to recompute using the new reference and 
perform a completely independent fit. Due to multiple scattering, the results 
of two fits using different reference zq's are not related by the track model, 
and cannot be obtained from one another. This effect is already present in a 
perfectly linear ideal case, where the detector resolution does not depend on 
x. It just reflects the fact that the weights of measurements to estimate the 
track parameters depend (eventually strongly) on zq. 

This is unfortunate, since in practice, track extrapolation is often desired from 
both ends of the track. Furthermore, optimal track estimates at every possible 
sensor position can be quite useful, either to detect outliers efficiently, or to 
optimally collect hits left over during a first pass. 

Finally, if one attempts such a traditional track fit in an experiment with large 
iV (such as NOMAD, where N ranges up to 150) iV can be too big to make 
inversion of practical. 

Given these difficulties, in NOMAD the Kalman filter was implemented in- 
stead [9]. 



4 The Kalman Filter 



The Kalman filter is a least-squares stepwise parameter estimation technique. 
Originally developed in the early 60's to predict rocket trajectories from a 
set of their past positions, it can be used to handle multiple scattering while 
estimating track parameters. We try here to briefly shed light on the features of 
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the Kalman filter for track fitting and refer to the literature for more details [1- 
4]- 

The Kalman filter technique gives, mathematically speaking, exactly the same 
result as a standard least squares minimization. In the framework of track fit- 
ting, it essentially avoids big matrix inversion and provides almost for free an 
optimal estimate of track parameters at any location, allowing the detection 
of outlying measurements, extrapolation and interpolation into other subde- 
tectors. 

The set of parameters x is called the state vector in Kalman filtering. 

Starting from Eq. 18 we rewrite the system equation in a stepwise form, where 
the state vector at location k is obtained from its value at previous location^ 
k — 1. 

x fc = f fc (x fc _ 1 ) +w fc (20) 

We shall assume in the following that both and Sk (the measurement 
errors from Eq. 15) are independent random variables with mean and a 
finite covariance matrix. 

Linearizing the system in the vicinity of Xfc_i, one obtains: 

ffc(Xfe-i) = F k ■ x fc _! (21) 
/i fc (x fc ) = H k ■ x fc (22) 

for the track model and the measurement equation. 

We can now recall keywords used in the Kalman filter estimation technique: 

• Prediction is the estimation of the state vector at a "future" time, that is 
the estimation of the state vector at time or position (k + 1) using all the 
measurements up to and including m k . 

• Filtering is estimating the "present" state vector based upon all present 
and "past" measurements. For Forward filtering, this means estimating track 
parameters at k using measurements up to and including m^. For Backward 
filtering, this means estimating track parameters at k using the measure- 
ments m 7v down to m k . 

7 One has to assume that the measurements are ordered with respect to time to 
handle multiple scattering because the covariance matrix of measurement residuals 
depends on this order. Without multiple scattering, the ordering does not affect the 
filter result as in the case of a traditional fit. 
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• Filter. The algorithm which performs filtering is called a filter and is built 
incrementally: filtering mj to consists in filtering rri\ to m^-i, propa- 
gating the track from rrik-i to and including m k . A filter can proceed 
forward (k increases) or backward (k decreases). 

• Smoothing means using all the measurements to provide a track param- 
eter estimate at any position. The smoothed estimate can be obtained as 
a weighted mean of two filtered estimates: the first one using mi to m k 
(forward), the other using to m k+1 (backward) 

One can understand the basic idea of the Kalman filter in the following way. 
If there is an estimate of the state vector at time (location) it is extrap- 
olated to time tk by means of the system equation. The estimate at time t k 
is then computed as the weighted mean of the predicted state vector and of 
the actual measurement at time tk, according to the measurement equation. 
The information contained in this estimate can be passed back to all previous 
estimates by means of a second filter running backwards or by the smoother. 

The main formulas for our linear dynamic system are the following: 

System equation: 

x fc = F k ■ x fc _i + w fc (23) 
E{w k } = 0, cov{w fc } = Q k (l<k<N) (24) 

Measurement equation: 

m k = H k -x k + £k (25) 
E{ek] = 0, cov{ £fc } = V k = Gk 1 (l<k< N) (26) 

where the matrices Qk and V k represent the process noise (multiple scatter- 
ing, bremsstrahlung, etc.) and measurement noise (detector resolution) re- 
spectively. The details of Q k calculation for the parameterization adopted in 
NOMAD can be found in Ref. [10]. 

As an example we include here the formulas for making a prediction: 

8 This leads to a subtlety in practice, when we actually have in hand the forward 
and backward filter estimates at k; averaging these would lead to double-counting 
the information from m k . Hence, to be proper, one must unfilter m k from one of the 
estimates. This small correction is implemented in our smoother and in the quantity 
Xk discussed below, but was deemed negligible and never implemented in the 
other breakpoint quantities. 
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• Extrapolation of the state vector: 

k—l T-i 

• Extrapolation of the covariance matrix: 
C k 1 = F k C k ^\F k + Q k 

• Predicted residuals: 
rJt 1 = m k - H^t 1 

• Covariance matrix of the predicted residuals: 

r>k— 1 t i I tt r~ik—\ ttT 

K k - V k + n k L, k H k 

Using the Kalman filter, the computer time consumed for a track fit is propor- 
tional to the number of hits on the track, while with the traditional technique 
it is proportional to the cube of the same number in case of multiple scattering. 

After the Kalman fitting procedure one has the following available informa- 
tion : 

• xf and xf : the Forward and Backward estimates of the state vector at 
position k, i.e., the estimate of the track parameters at location k using 
measurements 1 up to k (forward) or N down to k (backward). 

• Xk ^ an d Xk ^ : the minimum x 2 value of the forward and backward fits 
up to measurement k. 

• \/ <Xfc ' F - ) and V^ k ' B ^ : the covariance matrices of xf and xf respectively. 

• Xfc, Xtrk an d : the same quantities, determined from the smoothed 
estimates, the equivalent of a full fit done at location k. Note that the Xtrk 
minimum does not depends on the location at which the parameter are 
estimated (xlvk f° r ^fc i s independent of k). 

Thus, much information exists as the by-product of a track fit: at every hit 
on the track away from the ends, we have the results of three fits for the track 
parameters at that hit: a fit to the part of the track upstream, a fit to the part 
of the track downstream, and a fit to the whole track. This information is the 
input to breakpoint analyses. 



5 Earlier Applications to Breakpoint Searches 

A natural way to compare xf and xf is discussed by R. Friihwirth [3] and 
was implemented in NOMAD tracking [9,11] before developing our extensions. 
One simply constructs the x 2 °f the mismatch of all the forward-backward 
parameters at each hit k: 

Xl (FB) = - Xf f [V^B) + V ^)]-i (xf - xf ). (27) 
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The value of Xk ls easily computed from the following relationship which 
holds for any k: 



2 _ 2 (F) 2 (B) 2 (FB) /^x 

Xtrk — Xk + X k + Xk ■ I 28 J 



After the track fit, we find the hit k for which Xk is a maximum, i.e., 
for which the forward-backward mismatch in track parameters has greatest 
significance; we call this maximum value x 2 ^ FB ^ ■ One can assign a breakpoint 
at that k there if x 2 ^ FB ^ is above some threshold. Fruhwirth also investigated 
various combinations of x 2 ^ FB ^ with x\ Xk aR d degrees of freedom in 
the track fits, but concluded that x 2 < - FB * ) was his best breakpoint test statistic 
[3,12]. 



6 Some Additional Breakpoint Variables Based on Constrained 
Fits to the Specific Breakpoint Types 



One may suspect that previously defined breakpoint tests do not make optimal 
use of the available information. Any x 2 quantity is by definition insensitive to 
the arithmetic sign of differences, while in the processes of interest, a decrease 
in the momentum is expected. Furthermore, xt mixes all the parameter 
mismatch information together. The signed forward-backward mismatch in 
single quantities such as 1/R can be examined, but it has the problem that 
the other 5 parameters are not constrained to be the same. (Physical changes in 
1/R can result, for example, in fitted mismatches in as well as 1/R-) Finally, 
an optimal test should use a more fully developed breakpoint hypothesis, so 
that a more meaningful comparison of x 2 's, with and without breakpoints, 
can take place. 



6. 1 Constrained Fits to Breakpoints 



Here we show how to obtain and examine the result that one would get by 
doing a traditional fit which uses all hits, but which allows a subset of the track 
parameters to be discontinuous at a particular hit k. We parameterize the full 
track with 1 to 3 added parameters in order to incorporate breakpoints of 
Types I, II, and III at hit k. E.g., for Type I, we replace the parameter 1/R 
by two parameters, a forward value 1/Rf just before hit k and a back value 
1/Rb just after hit k. Thus, our fits to the 3 breakpoint types have 7, 8, and 9 
parameters, respectively. We denote these sets with breakpoints by a. rather 
than x, and they are, respectively for the three types: 
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cki = {x,y, 1/Rf, l/-RB,tanA,0,t}, 
an = {x, y,l/R, tan A^, tan \ B , <p F , <p B ,t}, 
"in = V, 1/ Rf, V Rb, tan A F , tan A B , (f> F , (f> B , t}. 



(29) 
(30) 
(31) 



For definiteness, we discuss here the concept in terms of a Type I breakpoint. 

One can imagine a cumbersome procedure whereby one puts a Type I break- 
point at a particular hit k, and performs a traditional x 2 track fit (with 7 
parameters in our case) to all the hits of the track, minimizing the full Type 
I track's x 2 , which we call 

xL,i,k(«i)- (32) 

One would obtain the best estimate of ai, its covariance matrix, and the min- 
imum value of Xfuiiifc) an f° r a breakpoint at that particular hit k. One could 
then repeat this for each possible value of k, obtaining numerous potential 
track-with-breakpoint fits. 

Essentially the same set of results can be obtained far more economically by 
starting from the results xf and xf at each k which already exist as by- 
products of the Kalman filter fit. These results carry all the information that 
we need, since their error matrices contain the information (up to linear ap- 
proximation) on how xt ^ an d x\ ^ change when the track parameters 
change. We need only perform a linear x 2 minimization in which {xf,xf} 
(now playing the role of the "measured data") is compared to the Type I 
breakpoint model prediction Hiol\. Hi is the model matrix with 12 rows and 
7 columns containing only zeros and ones : 

HjOli = {(x,y, l/-Rp,tanA,0,t), (x,y, l/R B , tan A, 0, t)}. (33) 

Or, introducing the two 6x7 submatrices Hf and Hf of H I : 

H iai = (H[ ai ,Hf ai ). (34) 

Since our 12 pieces of "measured data" {xf,xf } are two independent sets of 
6 parameters with their corresponding covariance matrices, the appropriate 
chisquare can be written as : 

Xl (FB) (") = («T - # F "f [V^Y 1 (Xf - H F a) 

+ (xf - H B a f [V^]- 1 (Xf - H*a). (35) 

(Here and below, we suppress the subscript I since the equations are true for 
all breakpoint types, with Eqns. 33 and 34 suitably changed.) 
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The full x 2 of Eq. 32 can be written as : 

xLi,i )fe («) = (m - h (H a )f [Vt m) ] 1 (m - h (//a)) (36) 

where vjf 1 ^ is the block diagonal matrix containing the covariance matrix 
y{m,F) q ^ ^ e measurements (m x . . . m k ) and the covariance matrix V k m ' B ^ of 
the measurement (m fe+1 . . . m N ). 

As shown in Appendix A, x\ ^ FB \a) of Eqn. 35 is related to xha.,i,k( a ) °f 
Eqn. 36 in a revealing way. At each hit k, 

2 / \ 2 (F) 2 (B) 2 (FBI/ \ , _x 



subject to sufficient linearity in the fits. Thus by finding the minimum of 

xl (<*)» we find the minimum of xLi,i,fc( a )> since xl (F) and Xfc (B) are 
known. 

The minimum of x\ ^ FB \&) is obtained without iteration since the model 
relating a to {xf, xf } is linear. The set of estimated parameters is given by : 

& k = v^H T (vf ) r 1 ^ (38) 



where X = {xf ,xf }; is the block diagonal matrix containing \/( Xfc ' F ) and 
y(x fc ,s). anc j w h ere the covariance matrix for the new estimate, V ( - (Xk \ is given 
by: 



yiPL k ) = 



H T (V k ^H 



-l 



(39) 



2 ( FB~) 

The value of at its minimum is 

xV FB \Z) = -(H T (vh- 1 x) T -*- (40) 



Since the elements of H are mostly 0, and the rest equal to 1, the multiplica- 
tions by H and H T were done by hand before coding the software; elements 
of (V k )~ l have only one or two simple terms. 

The computer time to perform the Type I, II, and III breakpoint fits at all hits 
k (away from the track ends) was a negligible addition (few per cent) to the 
NOMAD track finding and fitting software. This added time was more than 
paid back by the speedup in matrix inversion which was obtained by explicitly 
unrolling the loops in the DSINV routine from CERNLIB [13]. 
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6.2 Breakpoint Variables Based on the Constrained Fits 

From the wealth of information thus available at each hit, we discuss two 
of the most useful categories: 1) signed differences, in sigma, of breakpoint 
parameters, and 2) x 2 comparisons based on the Fisher F-test. 

As an example, for Type III breakpoint fits, we let Dm ;k (l/R) be the forward- 
backward difference in 1/R, divided by its standard deviation (sigma, com- 
puted from the covariance matrix taking account of errors in both quantities 
and their correlations). This signed quantity effectively gives the significance, 
in sigma, of the jump in momentum at that hit. Similar quantities, with analo- 
gous notation, are calculated for all components of a, for all breakpoint Types. 
Thus, for bremsstrahlung studies, Di k {l/R) gives the momentum change un- 
der the constraint that all other track parameters are continuous at hit k. 

The Fisher F statistic [14] is appropriate for testing if adding parameters yields 
a statistically significant reduction in the x 2 of a fit. It is simply the ratio of 
the respective x 2 /d°f for the two versions of the fits. Thus we naturally apply 
it to our track fits with and without breakpoint parameters. For example, for 
NOMAD's Type I fits, we have at each hit, 

Fi, k = (X 2 m /(N - 7)) / (xL/(N - 6)); (41) 
we similarly define F II fc and Fm^. 

For a true breakpoint at a given hit k, each F statistic in principle obeys the 
standard significance-table distributions [14]. However, we normally search 
through all hits in a track for the lowest value of F, denoted F. This lowest 
value does not of course follow the usual significance-table values. Hence, as a 
practical matter, we use data or Monte Carlo events to measure the distribu- 
tion of F and the effect of tests based on it. 



7 Effectiveness of the additional breakpoint variables 

The effectiveness of any breakpoint search is of course highly dependent on 
details of hardware (e.g., number and quality of the position measurements) 
and software (e.g., how many tracks with breakpoints are reconstructed as a 
single track). We present for illustration some experience with the NOMAD 
detector, using both simulated and real data events. 
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Fig. 2. Test of breakpoint search criteria using real data (muons producing 

2 ( F B) 

(5-electrons in the NOMAD setup). Comparison of breakpoint chisquare (Xh ) 
and normalized difference between curvatures in backward and forward directions 
((1/Rb — 1/Rf)/{gi/r)) for muons (solid line) and electrons (dashed line). In both 
electron distributions, the excess on the right is evidence of potential breakpoints. 

7.1 Electron identification and reconstruction 

Algorithms developed for electron identification and reconstruction have been 
checked under running conditions using 5-electrons produced by straight- 
through muons (5 GeV < p M < 50 GeV) crossing the NOMAD detector during 
slow-extracted beam between neutrino spills. This sample of selected electrons 
from real data can be used to check the subdetector responses compared to 
simulations and to tune breakpoint search criteria taking into account the 
effect of the drift chambers alignment quality (see Fig. 2). 

A special approach to deal with electrons emitting bremsstrahlung photons 
(Type I breakpoints) in the NOMAD detector has been developed [10]. If one 
has identified a reconstructed track in drift chambers as being an electron^], 

9 Transition Radiation Detector (TRD) is used for electron identification in the 
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Fig. 3. A reconstructed event from real data (the projection onto the yz plane, in 
which tracks bend) before an attempt to apply the breakpoint search algorithm. 
The track at the bottom was identified as an electron by TRD. The triangles are 
track extrapolations used to search for more hits and to match information from 
different subdetectors. 




r 



Fig. 4. The same event as before after applying a recursive breakpoint search al- 
gorithm. Two breakpoints are found along the electron trajectory and they are 
associated with two photons (dashed lines): one built out of a conversion inside the 
drift chambers fiducial volume and the other from a stand-alone cluster in electro- 
magnetic calorimeter. The bars on the right are proportional to energy deposition 
in the electromagnetic calorimeter. 

then hard bremsstrahlung photons can be looked for and neutral tracks can 
be created requiring further matching with preshower and electromagnetic 
calorimeter. A successful application of this approach to an event from real 
data is shown in Fig. 3 and Fig. 4. A recursive breakpoint search algorithm has 
been applied to a track identified as an electron. As a result two breakpoints 

NOMAD experiment. 
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have been found, each associated with an observed hard bremsstrahlung pho- 
ton (one of which converts to a reconstructed e + e~ pair). When applying the 
breakpoint search algorithm to an electron track one must keep in mind a 
potential problem related to a possible presence of several kinks along the tra- 
jectory (since it can bias the calculation of variables used for the breakpoint 
search). Details are in Ref. [10]. 

7. 2 Studies using simulated pion decays 

As another example of using the breakpoint variables described in Sec. 6.2, 
we compare several tests for detecting the Type III breakpoint (discontinuity 
of l/R, tan A, and </>) of the decay it — > fiu, and for locating the position of 
the decay. 

The results shown here are for Monte Carlo (MC) simulation [7] of muon neu- 
trino interactions in the NOMAD detector. From the collaboration's standard 
MC samples, we selected two samples of reconstructed tracks: 

• 25000 pions which did not decay and which were reconstructed. 

• 2000 pions which did decay via n — > fiv, and for which a single track was 
reconstructed consisting of hits left by both the pion and the muon. 

Neither of these two samples contains the pion decays which were broken 
into separate tracks by the track-finding algorithm, with a vertex assigned at 
the decay point. The selected tracks were chosen to have more than 20 hits 
(N > 20) and no backward looping, and with the MC pion decays within 
these hits. Track-finding mistakes (for example use of hits from other tracks 
in the event) were included, but to obtain the pion track, we required that at 
least 90% of the hits be correctly assigned. 

Figure 5 contains histograms of momentum x charge for the two samples. In 
order to have samples with similar momenta, we consider here only tracks 
with momentum less than 6 GeV. 

Figure 6 shows, for one of the decaying negative pion tracks, the values of 

2 ( FB} 

Xk j ^ni,fc) Dm t k(l / R) and Dm^{(p) at every hit where they are computed. 
They are plotted at the z-positions of the hits, which in the NOMAD detector 
are every few cm. The decay point in the MC is indicated by the line at 360 
cm. The dotted lines are at the ^-positions of the extreme values (maximum 
or minimum, as relevant for a breakpoint) of the respective variables. Near the 

2 ( FB} 

MC decay point, both Xk an d F\u,k reach their extreme values, indicating 
respectively: a large forward-backward mismatch in the Kalman fits, and a 
marked improvement in the track x 2 by adding a three-parameter breakpoint. 
The sign of the change in l/R corresponds to a decrease in momentum, as 
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Fig. 5. Histogram of momentum x charge for simulated pions without decay (left) 
and with decay (right). 

expected in a decay. Dm k {(f)) also shows an extremum near the breakpoint; 
however, we have generally not used it as the primary breakpoint indicator. 



7. 3 Tests for Existence of a Breakpoint 



For testing the existence of a breakpoint, we compare the following test statis- 
tics: 

• x 2 /dof for the breakpointless fit. (This test gives no additional information 
on the location of the breakpoint.) 

• X 2 ( FB ^. The breakpoint is located at the maximum value of the forward- 
backward mismatch chisquare (xk < " FB ' ) ) for breakpointless fit among hits in 
the track; 

• F. The breakpoint is located at the minimum value Fisher F among hits in 
the track (applicable separately to the different breakpoint types if desired). 

• Fni combined with the forward-backward mismatch in radius of the curva- 
ture Dm k {l/R) or angle Dm k {(f)). F m gives the breakpoint location and 
-Din,fc(l/i?)> -Dni,fc(0) are computed at that location. 

Fig. 7 contains histograms of x t 2 r k/ dof > Xf u ii,i/ dof > xLi,n/ dof > and xLi,ni/ dof 
for decaying and non-decaying pions, with the means shown. The x 2 / d °f is 
in general reduced by adding breakpoint parameters. A measure of the signif- 
icance of the improvement is the Fisher F (Eq. 41), which is shown in Fig. 8 
for the same fits. Also shown in Fig. 8 is a histogram of x 2 ^ FB \ the maximum 
forward-backward mismatch from the breakpointless fit. 
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Fig. 6. x 2 ( FB \ Fui t k, Dm t k(l/R) and Dui t k(4>) as a function of the z of the DC 
hit. The MC decay point is shown in solid line. 

From the histograms in Fig. 8, one can calculate the efficiency of labeling tracks 
as pion decays using "cuts" on these variables. For tabulating a comparison, 
we choose a cut value for F such that 10% of non-decaying pions are falsely 
called pion decays. We then find the efficiency for identifying real pion decays, 
i.e., the percentage of pion decay tracks which are called pion decays when 
using the same cut value. This cut value in principle depends on the dof, but 
for illustration we use the same cut value for all track lengths. The results 
are given in the first five columns of Table 1. We include for comparison the 
result for testing for pion decay simply by using Xtrk/d°f; i- e -> the test one 
would naturally use if one had only a traditional non-Kalman track fit. As 
observed by Fruhwirth [3], this test is quite competitive with the x 2 ^^-test 
for detecting the existence of a breakpoint, even though it gives no information 
about the location. The highest efficiency is obtained using F in - We find this 
to be true for various other comparisons, although one is cautioned that the 
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Fig. 7. Histograms of x? rk /dof, xf ullI /dof, Xfuii,n/ dof > and Xfuii,ni/ dof for the sample 
of pion decays (shaded) and pions with no decay (white) . The samples are normal- 
ized to the same number of events. 

particular efficiencies listed are for the sample of decaying pions which were 
not detected by the standard track finding/fitting, and are hence expected to 
be highly experiment-dependent. 

Next we add the signed information available in the new fits: the difference 
in radius of curvature Dm^(l/R), and/or the angular differences -Dni,fc(</>) 
and _D in ,fc(tan A). We recommend that Dm^{l/R) not be used to locate the 
breakpoint, but rather we evaluate this difference at the location dictated 
by Fm. Fig. 9 contains scatter plots of F m vs Dm tk (l/R) for decays and 
non-decays. Because true decays have a decrease in momentum, a judicious 
cut on this scatter plot is more effective than a cut solely on Fm- An even 
more powerful technique is to construct likelihood functions based on these 2D 
densities, and use the constructed likelihood ratio as a test of pion decay. The 
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Fig. 8. Histograms of the x 2 {FB) and F h F u , F m for pion decays (shaded) and 
pions with no decay (white). 



Table 1 

Efficiency of the detection of pion decays, using the various test statistics. The first 
five columns are for simple cuts on the respective variables. The last column is for 
a cut on a 2D likelihood-ratio test for Fm vs Dm t k(l/R) described in the text. The 
cut value for each test is chosen so that 10% of non-decaying pions are wrongly 
called decays. The efficiency is computed with respect to the sample which contains 
only tracks which were not detected as decaying by the original track-finding/fitting 
algorithms. 



X t 2 rk /dof X 2(FB) 


Fj Fu 


-Fin 


F m vs Dm,k(VR) 


45% 41% 


38% 40% 


49% 


56% 
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Fig. 9. Fni vs Dm^(l/R) and Fm vs Dm^(4>) histograms for the pion decays and 
pion with no decays simulated tracks. In order to superimpose negative and positive 
tracks on the same plots, we have changed the sign of Din^{l/R) for one charge. 
(Recall that 1/R is a signed quantity reflecting the charge of the track.) 

result of such a procedure (applied with simple smoothing of the 2D densities, 
and tested on an independent sample) is given in the last column of Table 1. 
The improvement over previous breakpoint tests [3] is most significant. Nearly 
as good efficiencies are obtained from scatter plots of Fm vs Dm t k((j)) (also 
shown in Fig. 9), and from F m vs Dm ^ (tan A). 



7.4 Finding the Location of the Breakpoint 



One may also ask which of the variables gives the best determination of the 
location of the breakpoint. We studied in particular the difference between the 
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Fig. 10. x 2 {FB "> and F m resolutions for MC pion decay sample in z. The white 
histograms are for the initial sample; the shaded histograms are for events in which 
detected breakpoints passed the selection criteria. 



MC decay point and the z position of the hit corresponding to the extrema 
X 2 or Fni- These are histogrammed in Fig. 10 for pions which decay. The 
white histograms show the resolution for the initial samples, with no selection 
made using x 2 ^ FB ^ or (-^ni vs Dm^{l / E)) . The shaded histograms show the 
(more relevant) resolutions for tracks remaining after decay selection using the 
respective variables. There was not a significant difference, within the limited 
scope of this study. 



8 Conclusions 



Replacement of mismatch chisquare for all the forward-backward parameters 
by the breakpoint variables introduced in Sec. 6.2 can give added power to 
breakpoint detection in the framework of Kalman filtering technique. We show 
in particular above that this is the case in a realistic simulation of pion decays 
in the NOMAD detector. In addition, these breakpoint variables have been 
successfully used to reconstruct electron hard bremsstrahlung in real data. As 
expected on theoretical grounds, our most powerful breakpoint detection is 
based on a scatter plot of a Fisher F-test vs. an appropriate signed difference 
of a track parameter across the breakpoint. 
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A Derivation of Eq. 37: Xf u \\ i 



2 (F) 

= Xk 



2 (B) 



+ Xk +Xk 



2 (FB) 



In Eqn. 36, vjf 1 ^ is the block diagonal matrix containing the covariance matrix 
vjf a,F ' > of the first k measurements m F = (mi . . . and the covariance 
matrix vjf 1 '^ of the last N — k measurements m B = (m fc+1 . . . m N ). The 
right-hand side of Eqn. 36 can thus be split into two terms: 



Xfull,I,fe( a ) — 



+ 



m 



m 



in' - h(H F a) 
in' 1 - h(H B a) 



(A.l) 



where one recognizes, in analogy with Eqn. 17, the forward and backward \ 2 
terms. We can expand each around their respective minima x. F and 5t B , and 
recall that covariance matrices are the inverse of curvature matrices, giving : 



*W«) = Xfc (F) + (Ax F ) T [\/^ fc ' F) ] _1 Ax F 
+X 2 k {B) + (Ax B ) T [V^]- 1 Ax B , 

where Ax F = xf - H F ct and Ax B = xf - H B a. 
Combining this with Eqn. 35 yields Eq. 37. 



(A.2) 
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